Optimal control of atom transport for quantum gates in optical lattices 
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By means of optimal control techniques we model and optimize the manipulation of the external 
quantum state (center-of-mass motion) of atoms trapped in adjustable optical potentials. We con- 
sider in detail the cases of both non interacting and interacting atoms moving between neighboring 
sites in a lattice of a double- well optical potentials. Such a lattice can perform interaction-mediated 
entanglement of atom pairs and can realize two-qubit quantum gates. The optimized control se- 
quences for the optical potential allow transport faster and with significantly larger fidelity than is 
possible with processes based on adiabatic transport. 

PACS numbers: 03.67.-a, 34.50.-s, 



I. INTRODUCTION 

Quantum degenerate gases, such as Bose-Einstein con- 
densates (BECs) [1] or cold Fermi gases [2|, trapped in 
optical lattices, provide a flexible platform for investigat- 
ing condensed matter models and quantum phase tran- 
sitions [3[. It has been proposed to use these systems 
as quantum simulators of solid state systems Q and 
for implementing quantum information processing (QIP) 
[E IE • Experiments on neutral atoms have shown some 
of the ingredients needed for QIP: the preparation of a 
Mott insulator state with just one particle per well, which 
is used as the initial state of a quantum register 0] , single- 
qubit rotation [8], and controlled motion of atoms so as 
to effect entangling interactions [8|, [9[ . 

A general requirement of QIP is accurate control of a 
quantum system. Often this includes control of degrees of 
freedom other than the qubit or computational basis, for 
example the center of mass motion of an ion or atom for 
which the spin (internal state) represents the qubit. One 
approach to achieving such accurate control is adiabatic 
manipulation of the relevant Hamiltonian. Unfortunately 
adiabaticity limits the speed of operations. One way to 
overcome this difficulty is to use optimal control methods 
0, E3] • Here we show that such techniques could improve 
the speed and fidelity of transport of atoms in an optical 
lattice. 

Recent experiments [E 0j], 03 have shown that quan- 
tum gates could be implemented in controllable opti- 
cal potentials by adjusting the overlap between atoms 
trapped in neighboring sites of an optical lattice. High 
fidelity of this dynamic process could be achieved by adia- 
batically changing the trapping potential. This, however, 
limits the overall gate speed to be much lower than the 
trapping frequency [7|, LL3[ . Here we present a detailed nu- 
merical analysis of the transport process used to effect a 
two-qubit quantum gate in [12[ , which is performed with 



the controllable double-well optical potential described 
in [lU and find that it gives an accurate description of 
the evolution measured in the experiment. Then we ap- 
ply optimal control theory to the transport process of 
the atoms, both with and without interactions, to show 
how to increase the speed of the gate. The success of 
this method in this specific problem demonstrates the 
promise of optimal control for coherent manipulation of 
a diverse class of quantum systems. 

II. THE EXPERIMENT 

A two-qubit quantum gate with neutral atoms can 
be realized in optical lattices through a controlled 
interaction-induced evolution of the wavefunction that 
depends on the states of the two atoms [E @ . Because 
atoms in their electronic ground states generally have 
short-range interactions, in order to use these contact 
interactions to produce entanglement, the atomic wave- 
functions must be made to overlap. Once the interaction 
has taken place for a fixed time, the two atoms can be 
separated again thus finishing the gate. In this paper we 
consider the control of such motion in a specific setup; 
however our theory can be applied to more general sys- 
tems. 



A. The Double- Well Lattice 

Neutral 87 Rb atoms are loaded into the sites of a 3D 
optical lattice obtained by superimposing a 2D optical 
lattice of double-wells [14| in the horizontal plane and 
an independent ID optical lattice in the vertical direc- 
tion. The horizontal lattice has a unit cell that can be 
dynamically transformed between single- well and double- 
well configurations. The horizontal potential experienced 
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by the atoms is [151 ] : 
V(x,y) = -V 



cos 2 I ^ ) (cos 2 ky + cos 2 kx) + 



sin 2 ( ^ ) (cos % + cos(kx — 0)) 2 



(i) 



where x and ?/ are the spatial coordinates, k = 2tt/X is 
the laser wave-vector and A is the laser wavelength. The 
potential (pp) depends on three parameters: 

(i) the strength Vo of the potential wells; 

(ii) the ratio tan of vertical to horizontal electric 
field components; 

(iii) the phase shift 6 between vertical and horizontal 
light components. 

The angle (3 determines the height of the barrier between 
adjacent double- well sites: by changing (3/ir from to 0.5 
the potential changes from a symmetric double- well con- 
figuration, with a spacing of lambda/2 (lambda/2 lat- 
tice), to a single- well configuration, with a spacing of 
lambda (lambda lattice). By changing (3 and together 
one varies the energy offset (tilt) of a well with respect 
to the neighboring one. The tilt of the double well is 
zero for 0/tt = ±0.5, while it is maximum (with a value 
depending on (3) for 0/ty = or ± 1. 

To effect a quantum gate, one varies the three parame- 
ters in time so as to move atoms occupying adjacent wells 
into the same well, allowing them to interact and finally 
returning them to their original positions. 

In Fig. [1] we show four snapshots of the cross-section 
of the optical potential along the direction of the double 
wells (x), and of the single-particle wave functions of the 
two atoms during a particular transport sequence. In the 
initial configuration each atom is prepared in the ground 
state of separate wells so that the properly symmetrized 
initial state is: 



|Vl} 2 ) (2) 



where ipL and ipR are wavefunctions localized in the left 
and right well, respectively, and 1 and 2 are the labels 
of the two (indistinguishable) atoms (see Fig. []Ji). For 
quantum gate operation, we would also include the inter- 
nal state of the atoms, but here we concentrate only on 
the external state (center-of-mass motion). ipL and i/jr 
are linear combinations of the lowest symmetric and an- 
tisymmetric energy eigenfunctions of the single-particle 
potential. In Eq. ([2]) and throughout the text we use the 
convention that single- (two-) particle states are labeled 
with lower (upper) case greek letters. 

The potential is changed by lowering the barrier and at 
the same time lowering the right well with respect to the 
left one, Fig. [D))-c). The atom initially in the right well 
remains in the ground state, evolving into the lowest state 
0o of the final potential, while the atom initially in the 



a) 

U -100 



-200 



1 1 


i 1 i 1 i 


1 i 










-0.4 -0.2 0.2 0.4 
xA + 0.25 



-0.4 -0.2 0.2 0.4 
xA+0.25 




-200 




-200 



-0.4 -0.2 0.2 0.4 
xfk + 0.25 



-0.4 -0.2 0.2 0.4 
xA + 0.25 



FIG. 1: (Color online) a) Initial configuration: the potential 
(solid lines) is in the A/2 configuration and the single-particle 
wave functions are localized in the left (dashed) or right (dot- 
ted) well. b)-c) Intermediate snapshots obtained by lowering 
the central barrier and tilting the potential, d) End of the 
process: the particle initially in the right well ends in the 
ground state of the single well, while the particle initially in 
the left well ends in the first excited state. 



left well evolves into the first excited state 0i, Fig. [TJi). 
When the two atoms are in the same potential well they 
interact through the usual contact interaction, which can 
be used to generate the entangling operation needed to 
realize a two-qubit quantum gate [j 



B. Experimental Procedure 

The experimental characterization of the transport 
process is accomplished by performing the potential 
transformation depicted in Fig. [T] with atomic samples 
loaded either only in the left sites or only in the right 
sites of the double wells [15[ . 

Briefly, Bose-Einstein condensates of 87 Rb atoms with 
4 • 10 3 < TVbec < 2 • 10 4 are loaded in the sites of 
the A-lattice with an exponential ramp of 200 ms du- 
ration. This loading populates only the ground band 
of the optical potential with mostly one atom per lat- 
tice site [l6|. Then the potential is transformed to the 
A/2-lattice in such a way that the atoms eventually oc- 
cupy either only the right sites or only the left sites of 
the double wells @, [l5[. Starting from this initialized 
state, we perform the transport process illustrated in Fig. 
QJi)-d). At the end of the process we measure the occu- 
pation of the lattice bands. To this purpose, we map 
the quasi-momentum of atoms occupying different vibra- 
tional levels of the optical potential onto real momenta 
lying within different Brillouin zones [I?], EH- This is 
achieved by switching off the optical potential in 500 /as 
and acquiring an absorption image of the sample after a 
13 ms time-of-flight. In this way atoms occupying differ- 
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ent vibrational levels appear spatially separated, allowing 
us to measure the amount of population in each vibra- 
tional state. 

The comparison between these measurements and the 
theoretical model requires an accurate determination of 
the evolution of the parameters Vb, /3 and characteriz- 
ing the optical lattice during the experimental sequences. 
The parameter Vb, which corresponds to the depth of 
the potential when it is set in the A/2 configuration, is 
measured by pulsing the A/2-lattice and observing the 
resulting momentum distribution in time of flight pjjj]. 
The parameters (3 and 0, which determine the shape of 
the double-well potential and are controlled using two 
electro-optic modulators (EOMs), are determined from 
measurements of the polarization of the laser beams af- 
ter the EOMs as a function of their respective control 
voltages [14j . 

We perform two series of experimental sequences in or- 
der to study the properties of the atomic transport as a 
function of the duration of the process and of the energy 
tilt between the two potential wells during the merge. In 
a first series of measurements the sequence involves con- 
verting the lattice from the double- well to the single- well 
configuration by changing /3, rotating the polarization of 
the input light using a linear ramp, while leaving constant 
the light intensity and the setting of the electro-optic 
modulator EOM# dedicated to the control of the phase 
shift 0. This sequence is repeated for different durations 
of the linear ramp from T = 0.01 ms to 1.01 ms. In a sec- 
ond series of measurements we consider the dependence 
of the transport on the tilt of the double-well potential 
during the merge. We perform the lattice transforma- 
tion using always the same duration of T = 0.5 ms, the 
same intensity of the light beam and the same ramp for 
changing the polarization angle /?, while the the setting 
of EOM# is kept constant in time during a sequence. We 
then repeat the sequence for different settings of EOM#. 
The time dependence for the three lattice parameters Vo , 
/3 and for measurements of the first series and of the 
second series are shown in Fig. [2^ and Fig. [2Jd, respec- 
tively. Fig. [2fc> shows the evolution of the parameter for 
two different settings of EOM#. The potential parame- 
ters are determined using our calibration of the lattice 
setup, taking into account effects such as different losses 
on the optical elements for different polarizations of the 
lattice beams and the dependence of the optical potential 
on the local polarization of the light [8] . These effects are 
responsible for the change of the potential depth Vo and 
of the angle 6 during the sequence despite the fact that 
both the intensity of the light and the settings of EOM# 
are not actively changed. 



III. THEORETICAL MODEL 

Here we describe the theoretical methods that we im- 
plement for investigating the dynamics in the system de- 
scribed above, starting with the case of non-interacting 




FIG. 2: Two possible sequences a (left) and b (right) employed 
to shift the atoms from a double- to a single-well configura- 
tion are shown in the left and right part of the panel. For 
each sequence we show the time-dependence of Vb, (3, for a 
sequence duration T. For sequence b we show the time de- 
pendence of for two settings of EOM#: —0.42 tt (solid) and 
-0.48 tt (dashed). 

particles. Then, we consider the experimental realization 
of the merging of adjacent lattice sites into a single site 
shown in Fig. [T] and we compare the results obtained in 
the experiment with the expectations based on our theo- 
retical model. This stage represents a useful benchmark 
to evaluate the reliability of the numerical model as well 
as for gaining insight into the details of the optical po- 
tential experienced by the atoms. Finally, we present the 
technique for optimizing the transport sequence, and we 
show how we can achieve a significantly higher fidelity at 
fixed operation time for the atomic motion than by using 
smooth sequences based on adiabatic evolution. 



A. Theoretical Framework 

We consider the ID problem restricted to the x axis 
by assuming that the optical potential can be separated 
along the three spatial directions, allowing us to express 
the atomic wavefunctions as a product of three indepen- 
dent terms. We consider the harmonic approximation of 
the potential in the y and z directions, having trap fre- 
quencies v y and v z respectively, that can be calculated 
as shown in [20[ and we assume that along y and z the 
atoms always occupy the lowest vibrational state. This 
restriction does not put limitations in studying dynamic 
processes involving low energy states of the double-well 
potential since it can be chosen to have non-degenerate 
vibration frequencies along all three directions, with the 
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lowest frequency always along x. We calculate the eigen- 
states of the system along the x direction by solving the 
eigenvalue equation using the finite difference method 
[211 ]. For the time evolution we consider the integration 
of the time dependent Schrodinger equation using the 
Crank-Nicolson method [22[. This method has the ad- 
vantage of being unconditionally stable and the error in 
the results scale quadratically with the number of space- 
time grid points in which the Schrodinger equation is 
solved. The relative error of the data presented is al- 
ways less than 10 -3 . In Appendix lAl we present a more 
detailed description of our numerical methods. 



B. Comparison to experimental results 

In this section we present the theoretical analysis of the 
transport processes described above and we discuss the 
agreement between the model and the experimental mea- 
surements. We start by considering the time evolution of 
the Hamiltonian spectrum during the two sequences a 
and b shown in Fig. [U 
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FIG. 3: Instantaneous spectrum of the ID Hamiltonian for 
sequences a and b for EOM#: —0.42 7r. 

At time t = the spectrum is made of nearly degener- 
ate doublets of almost equally spaced pairs of harmonic 
oscillator states, while at time t = T the levels are similar 
to those of a single harmonic oscillator 1 . 



Fig. [3] shows the time evolution of the eigenstates of 
the single-particle Hamiltonian. The atoms initially pre- 
pared in the two local ground states in the right and left 
wells (iJjr and i/jl) evolve into the instantaneous eigen- 
states ending in the ground and first excited state of the 
final configuration, respectively. This approach requires 
the sequence to be performed slowly with respect to the 
timescale associated with the gaps between the relevant 
energy levels. The optimal "speed" in the parameter 
space can be calculated using the Landau-Zener theory 
for avoided level crossings. 

For gaining quantitative insight into the properties of 
the transport we perform numerical simulations for the 
sequences used in the experiments, also taking into ac- 
count possible deviations of the parameters from the ex- 
perimental calibrations, and we compare the results with 
the experimental measurements. The relevant quantities 
for our analysis will be the overlap f% of the energy eigen- 
states (j) n of the final potential with the evolved state ip a 
where a = L,R indicates the initial well occupation: 



fz = \{<t> n \u{T) ivur 



(3) 



where the operator U(T) is the single-particle time- 
evolution operator from time t = to time t = T. In 
the experiment can be measured as the population of 
each energy level at the end of the process. 

We now consider how the atoms evolve when the pa- 
rameters change according to sequence a of Fig. [2] as a 
function of the total time T, focusing on atoms start- 
ing in |^l) 2 - In Fig. [H we show the final population of 
the ground /q', first and second excited state f£ mea- 
sured in the experiments and calculated for four values of 
which differ from the one of Fig. [2] a by a constant offset 
AO a 3 . The results obtained by the model are in reason- 
able agreement with the experimental observations; we 
find best quantitative matching for A6 a /7r = -0.03, for 
which the rms deviation between model and theory is 
reduced from 0.13(at AO = 0) to 0.08. 

Now we consider the sequence b of Fig. [2j performed 
for different ending values #5 of the parameter around 
— 0.477T. As shown in Fig.[5j both the experiment and the 
model show a strong dependence on #5 for the transport 



1 We have verified that the results for the ID spectrum are in good 
agreement with full calculations in 2D (restricted to states with 



vibrational excitation along the x direction). Additional energy 
levels are present in the 2D spectrum, associated with states 
with vibrational excitation along y. However, those states can 
be neglected for studying the dynamical process considered here 
since their energy is always higher than the three lowest states 
of the ID spectrum. 

2 Both in the experiments and in the simulations the evolution of 
the atom initially in the right well, i.e. in state \i/jr), shows a 
weaker dependence on the properties of the sequence and is less 
instructive. For instance, in the simulations for T = 0.5 ms the 
population in the ground state is of order of 99% for a broad 
range of parameters. 

3 We do not consider variations in Vb and (3 due to the small 
dependence of the transport process on those parameters within 
the range associated with the accuracy of our calibrations. 
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FIG. 4: Population of the first three eigenstates of the Hamil- 
tonian, ground (top), first excited (middle), second excited 
(bottom), at the end of sequence a as a function of the se- 
quence duration T. The experimental data (symbols) are 
compared to the four sets of numerical data (lines) obtained 
for 0/tt = -0.454 + A^/tt, with A^/tt = (dot-dash), - 
0.02 (dash), -0.03 (solid) and -0.04 (dot), while -0.454 is the 
nominal value of 0/tt expected from the calibrations. 



FIG. 5: Population (overlap absolute squared) of the first 
three eigenstates of the Hamiltonian at the end of sequence 
b as a function of 0b- The duration of the sequence is fixed 
to T = 0.5 ms. The experimental data (symbols) are in good 
agreement with the numerical data (lines). In this graph the 
x axis for the experimental has been shifted by an offset of 
-0.015 with respect to the initial calibration. 



ibration by an offset AO^/ir = -0.015, which reduces the 
rms deviation from 0.4 to 0.15. The axis for the exper- 
imental data in Fig. [5] has been corrected by the offset 
AOfr. Thus, while showing the reliability of the model in 
describing the dynamic process, the comparison between 
theoretical and experimental results also allows one to re- 
fine the calibration of the parameters characterizing the 
optical potential. 

Finally we find that adding an offset of AO/tt — -0.016 
to the calibration brings the data from both sequences 
to a good agreement with the theory and reduces the 
rms deviation from 0.19 to 0.11. This is three times 
larger than the uncertainty of the offset from our EOM 
calibration but is still consistent with measurements of 
the lattice tilt from [16j]. This discrepancy might be ex- 
plained by the birefringence in the vacuum cell windows, 
which is not accounted for in our model. Inclusion of this 
offset should improve both the predictivity of the model 
and the experimental optimization of the collisional gate 
based on the numerical technique described below. 



IV. OPTIMIZED TRANSPORT 



of the atom starting in the left site of the double well. 
The transport into the first excited state has a maximum 
theoretical fidelity of 0.95 for O^/ir = -0.474. Less nega- 
tive values of #5, i.e. increasing tilts, lead to a decrease of 
fidelity due to the increase in the fraction of population 
ending in the second excited state. Values of 6b /tt closer 
to -0.5, i.e. more symmetric configurations of the double 
well, lead to decrease of fidelity associated with larger 
fractions of population ending in the ground state. Also 
in this case the experimental data and the theoretical 
model are in satisfactory agreement. For these data the 
deviation between theory and experiment is more sen- 
sitive to the value of the phase shift #5. We find best 
agreement by shifting the value determined from the cal- 



In this section we employ optimal control theory to 
obtain fast and high-fidelity gates. Our aim is to find 
a temporal dependence of the control parameters Vfo(t), 
/?(£), 0(t) that improves the fidelity even for a shorter se- 
quence duration, when the adiabatic sequences presented 
above yield a poor fidelity. Quantum optimal control 
techniques have been successfully employed in a variety 
of fields: molecular dynamics [k| [23l 24], dynamics of 
ultracold atoms in optical lattices [25|, [26|, |27[, imple- 
mentation of quantum gates @, l2q|. 

We use the Krotov algorithm [29| as the optimization 
procedure. The objective is to find the optimal shapes of 
the control parameter sequences that maximize the over- 
lap (fidelity) between the evolved initial wave function 
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and a target wave function. The initial and target wave 
functions are fixed a priori. The algorithm works also 
for more than one particle. The method consists in it- 
eratively modifying the shape of the control parameters 
according to a "steepest descent method" in the space of 
functions (for more details see [7]). The method requires 
evolving each particle's wave function and an auxiliary 
wave function backward and forward in time according to 
the Schrodinger equations. In our simulations we use the 
Crank-Nicolson scheme to realize this step as described 
in Appendix lAl 



A. Non-interacting case 

We optimize the gate for T = 0.15 ms 4 choosing as a 
starting point for the optimization a sequence similar to 
Fig. [2]d where is for simplicity taken constant to the 
final value Ob/ir = —0.474. Without optimization the fi- 
delities for the atom initially in the left and right well are 
fi — 0.57 and f$ = 0.69, respectively. The infidelities 
are shown in Fig. [6] as a function of the number of opti- 
mization steps: the algorithm of optimization is proven 
to yield a monotonic increase in fidelity [l0[, however it 
does not guarantee to reach its 100% value. The results 
for the two atoms give a fidelity above 98.7%. 
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FIG. 6: Infidelities (1 — /") for the atom initially in the left 
(a = L,n = 1, squares) and in the right well (a = R, n = 0, 
plus) as a function of the optimization step. 

The resulting optimized parameter sequences are 
shown in Fig. [71 and compared to the original sequence 
without optimization. We find that the optimized se- 
quence for the potential depth Vq differs negligibly from 



4 We chose this time in order to show the benefits of the optimiza- 
tion procedure for a sequence duration which cannot provide a 
good fidelity with smooth parameter ramps based on adiabatic 
evolution. While a shorter duration could be chosen in principle, 
this choice can be easily experimentally implemented with no 
major changes in the present experimental apparatus, allowing 
future prosecution of our studies on this subject. Increasing the 
total time should improve the best fidelity. 
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FIG. 7: Initial (dotted) and optimized waveforms (solid) (3(t) 
and 6(t) as a function of time for a sequence of T — 0.15 ms. 



the initial guess. In principle, the algorithm could achieve 
still higher single-particle fidelities from different starting 
points. 

In Fig. [8] we show the square absolute value of the wave 
functions of the two atoms as a function of time, the 
ID potential time dependence and the projections of the 
initially left- well state onto the lowest four instantaneous 
energy eigenstates \<j) n (t))\ 



P n (t) = \{<Pn{t)\U{t)\^ L )Y 



(4) 



Notice that p n (T) = f%. As can be easily seen, the 
optimal time evolution is much less smooth than the adi- 
abatic one as it takes advantage of quantum interference 
between non-adiabatic excitation paths to obtain better 
results. 



B. Interaction effects 

Up to now we have considered only the single-particle 
evolution of the system, i.e. without including any in- 
teraction between the particles. This approximation is 
valid in our transport sequence as long as the two wave 
functions in nearby wells do not overlap. When the two 
particles overlap in the same well we must take into ac- 
count interactions, and we model them with an effective 
ID contact potential: 



Vint (|^i - x 2 \) = g\D&(xi - x 2 ) 



(5) 



where X{ are the coordinates of the two atoms and guj 
is an effective ID coupling strength [30] expressed by 
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FIG. 8: (Color online) Comparison between the evolution of the atoms with and without optimal control. Top (left to right): 
non optimized case, absolute square value of the wave functions as a function of time (atoms initially in the left and right 
well respectively); ID trapping potential as a function of time; projections p n (t) at time t of the state initially in the left well 
onto the instantaneous eigenstates \<p n (t)) with n — (blue solid), n — 1 (red dashed), n — 2 (green dotted), n — 3 (magenta 
dot-dashed). Bottom: analogous plots for the optimized case. 



9id — Zossh^iSyVzi where a s is the scattering length for 
87 Rb atoms and h is the Planck constant. The spectrum 
is modified by the interactions: the state with one atom 
in each well is lower by ~ 3 kHz than the doubly-occupied 
states. 

As in the case without interactions we start with each 
atom localized in a separate well. Notice that we are 
considering wave-functions that are symmetric under the 
exchange of the coordinates of the two particles. We 
consider the two-particle fidelity: 

2 



tg 



^int(T) 



(6) 



where U{ n t(T) is the two-particle evolution operator for 
the Hamiltonian of the two atoms, which includes inter- 
actions, ^in^ is an eigenstate of the two-particle Hamil- 
tonian at time t = 0, corresponding in the limit gw — » 
to the symmetrized product of the single-particle wave- 
functions localized in each well (see Eq.([2|)); the target 

state ^tg^ is an eigenstate of the two-particle Hamil- 
tonian at time t = T which, in the limit of vanishing 
interactions, corresponds to the state 

1 



l*tg> 



V2 



(l0o)il01> 2 + |0l)il0O> 2 ) 



(7) 



The square modulus of 



) and 



tg 



in the two- atom 



coordinate representation, are shown in Fig. [9^-b. In 



order to make a comparison between the interacting and 
non interacting cases we define a two-particle fidelity also 
in the non- interacting case: 

|2 



F=\($ tg \U 1 (T)®U 2 (T)\* i . 



(8) 



where U\{T) and ^(T) are the single-particle evolution 
operators for the two atoms without interactions. In Ta- 
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fi 


F 


i*int 


non optimized 


0.69 


0.57 


0.22 


0.22 


transport optimized 


0.99 


0.99 


0.98 


0.93 


interaction optimized 


0.98 


0.98 


0.96 


0.97 



TABLE I: Results of our numerical simulations for three dif- 
ferent sets of control parameters: the non optimized case 
Fig.[2]D; the transport optimized case Fig.[7]where the optimal 
control algorithm is used without taking into account interac- 
tions; the interaction optimized case where the optimal con- 
trol algorithm is used taking into account interactions. The 
quantities shown are: the single-particle populations and 
fx calculated without interactions, the two-particle fidelities 
F and Fmt calculated without and with interactions. 

bleUwe summarize our results for T = 0.15 ms obtained 
with three different sequences: first, the non optimized 
sequence Fig. [2Jd; second, the transport optimized case 
Fig. where we used the optimal control algorithm to op- 
timize the single-particle populations not taking into ac- 
count interactions; third, the interaction optimized case 
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FIG. 9: (Color online) Absolute square values of the relevant 
symmetric wave- functions in the coordinates of the two atoms: 
a Initial wave function in the state |^m^ with one atom in 
the left well and one atom in the right well; b wave function of 
the target state | i>t g ^ . c evolved wave-function using the non- 
optimized sequences of Fig. EJd giving a fidelity F int = 0.22 
(for T = 150/xs); d evolved wave- function using the optimized 
sequences of Fig. [3 giving a fidelity Fi n t = 0.93. 



where we apply the optimal control algorithm using as 
the initial guess the transport optimized sequence Fig. [71 
and then optimizing including the interactions in the evo- 
lution. 

The resulting wave- functions for the non optimized and 
transport optimized sequences are compared in Fig. Efc- 
d. Without optimal control the two-particle fidelity with 
and without interactions is F = F[ nt =0.22 while with 
(non-interacting) optimization we obtain F c± f± — 
0.98 and F[ nt = 0.93. This shows that interactions 
spoil slightly the efficiency of the transport process as 
one might expect. Optimal control can subsequently 
be applied while including interactions in the optimiza- 
tion, producing a control sequence with a fidelity of 

F int =0.97. 

Another consideration is the experimental bandwidth 
available for feedback control. The optimized control 
waveforms Fig. were obtained with no restriction on 
the frequency response of the control, and typically have 
frequency components on the order of a few times the 
lattice vibrational spacings (see Fig. [T0|) . i.e. larger than 
the bandwidth of our control electronics. Clearly, using 
a filtered version of these waveforms will lead to lower 
control fidelity and it will be important to increase the 
experimental bandwidth of the control electronics (cur- 
rently about 50 kHz). In addition, it may be useful to 
develop an optimization sequence that includes the lim- 
ited control bandwidth, although it is likely that frequen- 
cies on the order of the vibrational spacing will always 
be needed. 



T 1 1 1 I I I I 




2 3 456789 

1000 



Frequency (kHz) 

FIG. 10: The normalized Fourier transform magnitudes 
\/3(f)\ (solid) and \0(f)\ (dashed) of the optimized control 
sequences /3(t) and 0(t) shown in Fig. [7] The spectra are 
normalized to the value at the fundamental frequency 1/T — 
6.67 kHz. 



V. CONCLUSIONS 

We have presented a detailed, numerical analysis of the 
transport process of neutral atoms in a time dependent 
optical lattice. We show how to improve the fidelity of the 
transport process for T = 0.15 ms from F mt = 0.22, us- 
ing simple adiabatic switching, to F[ nt = 0.97, using op- 
timal control theory. We expect better results for longer 
control times. We analyze the effect of atom-atom inter- 
actions on the transport process and we show that the 
optimal control parameter sequences found in the non- 
interacting case still work when including interaction. We 
obtained the same transformation as in the case of the 
adiabatic transport with a better fidelity and in a time 
shorter by more than a factor of three, which represents a 
relevant improvement in terms of scalability of the num- 
ber of gates that can be performed before the system 
decoheres due to the coupling to its environment. This 
technique can be easily adapted to other similar trans- 
port processes and also extended to atoms in different 
magnetic states, which can allow the implementation of 
a fast, high-fidelity quantum gate in a real optical lat- 
tice setup with the qubits encoded in the atomic internal 
states [12|. In the future, it would be interesting to study 
the possibility of including the effect of errors in the op- 
timization procedure and thus investigate in more details 
the robustness and noise-resilience of the optimal control 
technique. 
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APPENDIX A: NUMERICAL METHOD 

In our numerical simulations we employ a finite differ- 
ence method (see for example [21]) that consists in dis- 
cretizing the coordinate representation in a homogeneous 
n points mesh in the interval [Xl;X2]: x k — x k -i = dx, 
xq = Xi . x n = X 2 . The number dx = (X 2 — X\ ) jn is 
the lattice spacing. In this discretized representation the 
eigenvalue equation becomes: 

(v(x k ,0) - e-^j ip a (x k ) = E^ a (x k ) (Al) 

where e = 3.5/(27r) 2 kHz, (3.5 is the conversion coeffi- 
cient between kHz and lattice recoils) and Vv(^fc) is the 
discretized wave function. The discretized second order 
derivative operator S 2 acts on any function as: 

S 2 x f(x k ) = f(x k+1 ) - 2f(x k ) + /(x fc _i). (A2) 

Expression (|A1|) is second order in dx 2 . If one arranges 
the function i/ja(x k ) in an n-dimensional array, then 
Eq. (|Aip can be rewritten as an eigenvalue problem with 
a n x n Hamiltonian matrix H. Since this matrix is tridi- 
agonal, the principal diagonal being V(x k ,0) — 2e/dx 2 
and the two sub-diagonals being filled with e/dx 2 , one 
can take advantage of its sparse structure for storage 
and computing. We calculate low energy eigenstates by 
approximately diagonalizing the Hamiltonian using the 
Jacobi-Davidson method [31]. This method is an iter- 
ative method which is capable of finding a few (typi- 
cally less than 10) eigenstates close to a chosen target 
energy. The advantages of using this and similar algo- 
rithms (Lanczos, Arnoldi) instead of exact diagonaliza- 
tion are that the method is much faster and one does not 
need to store the whole Hamiltonian. In our eigenvalue 
problem we take full advantage of this method, given the 
sparse structure of the Hamiltonian H. The error of this 
approximate method compared to an exact diagonaliza- 
tion method is negligible for our purposes. 



To study the time evolution of the wave functions of 
the atoms we integrate numerically the time dependent 
Schrodinger equation. Introducing a time slicing in the 
interval [0; T] with time interval dt, the Schrodinger equa- 
tion has the form: 

ip(x k: t n+1 ) -i/j(x k ,t n ) = -i dt H{x k ,t n )ijj{x k ,t n ) (A3) 

The discretized expression (|A3|) gives an iterative rela- 
tion to compute the wave-function at time t n +i, from 
the expression of the wave-function at time t n . This 
is one example of an explicit method: the coefficients 
ip(x k ,t n +i) can be directly calculated from ip(x k ,t n ). 
Explicit schemes have the great advantages of being ex- 
tremely fast and easily implemented. However this ex- 
pansion is only first order in dt and is not always stable. 
Therefore, we used the Crank-Nicolson scheme [22|, an 
implicit method, which consists in taking a time average 
of the right-hand side of (|A3ft between time t n and £ n +i, 
namely 

idt 

i/j(x k ,t n+1 ) - ip(x k ,t n ) = -—[H(x k ,t n )i/j(x k ,t n )+ 

+ H(x k ,t n+1 )^(x k ,t n+1 )] 

(A4) 

This method is of the second order in time and space 
and it is unconditionally stable. The price for all these 
advantages is that a tridiagonal set of linear equations 
must be solved to get V^n+i) as shown in Eq. (|A4|) . 
We used common Fortran routines to solve the linear 
equations problem [32^ 

We solve a 2D time dependent Schrodinger equation 
in the two coordinates of the atoms by making use of 
the extension in two dimensions of the Crank-Nicolson 
method called the Peaceman-Rachford method [21]. This 
is an implicit method and the integration proceeds in two 
steps: first the initial wave- function is integrated in time 
considering only one direction in the coordinate space, 
then from the intermediate wave-function we obtain the 
final wave-function by integrating in the other direction. 
This method is an example of alternate direction implicit 
schemes. 

In our simulations we used n = (X 2 — X\)/dx — 10 3 
and riT = T/dt — 5 • 10 3 that assures convergence of the 
results with a relative error which is less than 10 -3 . 
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